{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# One-way ANOVA\n",
    "\n",
    "For the one-way ANOVA, you can use either the function in \"scipy.stats\". Alternatively, you can use the \"statsmodel\" tools: they provide more, and nicer formatted, information.\n",
    "\n",
    "Author:  Thomas Haslwanter, Date:    Feb-2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "import urllib\n",
    "from statsmodels.formula.api import ols\n",
    "from statsmodels.stats.anova import anova_lm\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# The importance of the variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFJCAYAAACsBZWNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wVWV+B/DvTS5v4UZDhgtlNwnGEBI21kZ8250GsIs0\nC9h1iVCMVDoDS1c7Tou7RYRWpKvystSt467rrnY7TjN2OgtmcFUgte4gpVvYNBitSEKMqwlsJ0QJ\na07IhoSc/hHu5dybe+55e8779/OPkpuc89yHcH/n93veIrIsyyAiIiLH5bjdACIiorBiECYiInIJ\ngzAREZFLGISJiIhcwiBMRETkEgZhIiIil0SdvmFvb7/Q602bloe+votCr+ln7I9U7I9U7I9U7I9U\n7I9UIvsjHs/P+HXfZ8LRaK7bTfAU9kcq9kcq9kcq9kcq9kcqJ/rD90GYiIjIrxiEiYiIXMIgTERE\n5BIGYSIiIpcwCBMREblE1xKlFStWIBaLAQCKioqwc+fO5Gtvv/02nnvuOciyjKqqKjz++OOIRCL2\ntJaIiChANIPw0NAQZFlGQ0PDuNckScKePXvwL//yLygsLMSLL76Ivr4+FBYW2tJYIiKiINEsR7e1\ntWFwcBDr1q3D2rVr0dramnztnXfewdy5c7F7927cd999mD59OgMwERGRTpqZ8OTJk7F+/XqsWrUK\nH3/8MTZs2IBDhw4hGo2ir68Px48fx/79+5GXl4c1a9aguroapaWlqtebNi1P+AJotZ1Iwor9kYr9\nkYr9kYr9kYr9kcru/tAMwqWlpZg9ezYikQhKS0tRUFCA3t5ezJo1CwUFBfj93/99xONxAMAtt9yC\nU6dOZQ3CordEi8fzhW+F6Wfsj1SB7w9JQrT9FEYq5gFX5m1kE/j+MIj9kYr9kUpkf5jetnLfvn3Y\ntWsXAKCnpweSJCWDblVVFU6fPo3z589jZGQE7777LubMmSOkwUSkQZIwrfYOTFu6GNNq7wAkye0W\nEZFBmpnwypUrsWXLFtTX1yMSiWDHjh1oaGhASUkJFi9ejO985zv45je/CQD42te+hrlz59reaCIC\nou2nEO04Pfb/HafHMuKbb3W5VURkhGYQnjhxIp5++umUr82fPz/5/8uXL8fy5cvFt4yIshqpmIeR\n8rmIdpzGSPncsZI0EfmK40cZEpEgsRj6mg4bGhMmIm9hECbys1iMJWgiH+O2lURERC5hECYiInIJ\ngzAREZFLGISJiIhcwiBMRETkEgZhorCQJOD4ce6sReQhDMJEYXBli0t8+cvc4pLIQxiEiUIg0xaX\nROQ+BmGiEEhscQmAW1wSeQh3zCIKgytbXMbPdaFvRgm3uCTyCAZhorCIxYDS2wGeF0vkGSxHExER\nuYRBmIiIyCUMwkRERC5hECYiInIJgzAREZFLGISJKDNJQrSlmbtrEdmIQZiIxruyzeW0pYsxbclC\nRI8eYTAmsgGDMBGNk7LNZeeHmFZ3F/ecJrIBgzARjaPc5jKBe04TiccgTETjXdnmsq/xdYyUzQGA\nsf8ODjqXDXNMmkKAQZiIMovFMFKzEH1vHkFf4+sA4FxZWjkmzTI4BRiDMBFlF4sBU6Yg2vkhAGfK\n0jx6kcKCQZiINDl9FCKPXqSw4ClKRKTtyhhxtP3UWEC0+yhEp+9H5BIGYSLSJxbDyM23Bvd+RC5g\nOZqIiMglDMJE5C9cukQBwiBMRNY5FRgzLV1iUCYfYxAmSuCHuTkOrukdt3Sp9QTXE5OvMQgTAdwc\nwgIn1/SmL11K3NOJexPZgUGYCNwcwgpH1/QmttM8+Bb6mg5jpHq+uXuz6kEewSVKRLgaSKIdp7k5\nhFEurCFOLl2SJPTv/j4AYKR6vr57X6l6RDtOA5WVwIFfcB0yuYZBmMJHksYHjLBtDpGpD6xwY02v\nIpiOlM9FX9NhXT+mrHqgrW2sH7gemVzCcjSFS7ax30QgCUEADsL4t9khhJRjGisrWfUgVzEIU6hw\n7Dc4fWB6LFoxrozm5uA/dJGnMQhTqFieRBSACT2BORwhbZKWoWAalqoHeR7HhClcrIz9ZhqD9OOH\neJDGv7XGokWPfRMJxkyYwsdkFhSUMi6AcGSCARn7pmBjECbSKTBl3JAI1EMTBZaucvSKFSsQu/LE\nXFRUhJ07d6a8Pjo6ir/4i7/A4sWLUV9fL76VRF4QpDJuEGiUmrn2m/xAMwgPDQ1BlmU0NDSofs8z\nzzyDzz//XGjDiDwj7cOea0o9QM/4PB+ayAc0y9FtbW0YHBzEunXrsHbtWrS2tqa8fujQIUQiESxY\nsMC2RhK5xo5xxQDMsHab7lJzGMa+ydcisizL2b6hvb0d7777LlatWoWPP/4YGzZswKFDhxCNRnH6\n9Gk8++yzePbZZ/Hcc89h+vTpmuXokZHLiEZzhb4JItscPw58+ctX/3zsGHD77eavJ0nArbcCbW1j\nWybauU5VkoCTJ4Gqqsz30Hrdy5zsRyIbaZajS0tLMXv2bEQiEZSWlqKgoAC9vb2YNWsW9u/fj56e\nHvz5n/85zp49iwkTJuCLX/wiFi5cqHq9vr6LQt9APJ6P3t5+odf0M/ZHKsv9MaME0xTjin0zSgAL\n14u2NGNaW9vYH9ra0Hf0V/rL20aW26iUa5P9YWW5lVeW/Rz4xdV2DMrAoPG/F/57ScX+SCWyP+Lx\n/Ixf1wzC+/btw+nTp7F9+3b09PRAkiTE43EAwCOPPJL8vh/84AeYPn161gBM5DuCxxVNTxYyGDQz\nlWuVwV7rdVHtsBXH5ykANMeEV65cif7+ftTX1+Phhx/Gjh070NDQgLfeesuJ9hG5T+S4onKXp8Y3\nxsYydYwNG11uo7WcyuxyKy77IRJLc0xYNNGlDpZPUrE/Unm2P4xmlGYy0Axl45T+MFNW9lImLIBn\nfz9cwv5I5YlyNBGJZ7gcbKYsrlWuNVPOdXPZj1fGookE4o5ZRC4wVQ62a7mN0SVTbiz70btUjMu/\nyGcYhIkA5z+8rZwAJJJP9lceVzloPTH+m3zyXoiUGISJ3Prw9sBGEn6ZaDVSMQ8jZXOSf87ftHHc\n35Nf3guREoMwBZckjW22oRFUM354u1nWdPDenjiUQs/7jcXQv+eZ5B+jnR+OC7KeeC9EBnFiFgXT\nlewWHacxTWMW77i1u0Ul7s0Adnr2sdv7Kxt4vyPV87OvsXb7vRCZwEyYAslQaTJtfDZ6piv7z9qY\nqbpSUnWxLG7l7yljez1Q4icygkGYAiljaTJb8FR8eGcta4oeP05rk/CSqsdnCxt+vwyyFDAsR1Mw\nXcma4ue6xvZ7BvSXebOUNU1v95iJSilWWEk10/VVNgxwDUvIFHLMhCm4YrGxE49iMeNlXpWMS2Sm\nqtomQdmeb2YLM7ulEGMQplAQFjwFru+1ezZvYGYLe7ykTmQFy9EUDiLLnqJO77G7FOunUq/alpQB\n26uaKB0zYQoPL5Y9rbRJ5/paz73ndFkmu/mmpE5kEoMwkR95fYtGAyXkbIE2MCV1IhUMwkQ+5OkM\n0eADQtZAa2UMnmPJ5AMMwkRWufBh7+UM0cxM9KyB1kxJ3euVAqIrGISJrHDx8AdPnMKUgReOafR0\npYBIgUGYyAJXP+yNBC4ns3UPPCB4uVJApMQgTGSBrR/2ogKnIlvHrbc6FohdnZWd/iAAcHyYPIlB\nmMgKu7I+ZZl7yUJEjx7RF0AyBG5lto62tvCUZhXruTk+TF7FIEyUjUtrcVPK3J0fYlrdXdoBRGV8\nWpmto7IydKVZjg+TlzEIE6lxcYZtSuC8QiuAZNuLOpGto7nZ/hKxx5YGcXyYvIxBmEiF25Ou+poO\no6/xdYyUzQGgHUC01tvaNkarDLpuLg1SC/4emChGpIZ7RxOpSAS1xL7FjmdQsRhGquejf88zY+2p\nnp89gLixV3Ta3s7925/KfNSj2t7QNrVjXLAVtd83kWAMwkRqnApqRg4v0NFmK+cbG32v6dWCax55\nOPnaSNmcsWs5cAiD0HOeiRzEcjSFg9lxSruX2njl8AKTZeSUEnhxMXLPnkm+lsjgJ73aaPv74Lgv\n+RWDMAWfWoDxwAQi2w4vMPjeTAd85XjrgV+ktre8AtNq78A1Dz8EecJEc+9DL477kk+xHE2Bl7FU\nWTHPE+fUZh13NlsON1H+tTT+rSiBK9ur7PfI8CV8vutpDP1pvX39zHFf8iFmwhR4mTJKz6wd1drZ\nyUQ53NR7E5VJKto7UjEvObMbAPJefN7cNYkCjEGYgi9DgPHUGKLgnZ1MvzfR49+xWHJcGBjbdIQb\nZRClYhCmcEgPMEYzP7UxVoHjysKyczvGR9Pfp873PVI93zsPO0QexCBM4aU388sysUtX5qo3YInM\nztXem5mHhvT32dOjP2PnhCmirBiEiTSoZai6MlcjS3/sDliSBNx6q+GHhvT3Oek/moxl7Nkedkxm\n2ERBwSBMpCFjhipJwOCg5paShkvMNq5LjrafAtrasrclw0ND+vsfurNWzNIpKxm2GQzw5EFcokSk\nJX2pEHB1CVDZnLH9ncsrMi4lGrf0p6gE0ZZm57aVVBipmAdUVgJtbeOD55XdsjA4mHHnqfSlUiKW\nTvXv/r5mhi1syZHy3sXF6DvwC2DmTDHXJrKAQZhID8UM5mhLc8oxgwAwrW75uO0lE0EqGbCKSsZ/\nn5OBOBYDmpvRd/RXqcFTGaDK5mCkbA6inR+mBur0Nbgm1uSmVwUApDygJDJsO/bqTrl3dzcKly3G\n+bePAfF8YfcgMoNBmHSRJKC9PQcVFaOhn1uTnt0CSM3gWk8gf/O3U4LtyM23Inr0iPv7G2cInuln\nF/c1vg5MmSI8Wx9XFaieLybD1nvv4mJEu7sBALndXWOZf+ksYfcgMoNBmDRJElBbm4eOjlyUl19G\nU9PFcAfiDOXpbEF50quNGLqzFvmbNiYvkTzcwKjEIQtFJYie6RISrDIFR13XNHrgg0qQtZph6xKL\noe/AL1C4bDFyu7u4XIo8g0GYNLW356CjIxcA0NGRi/b2HNx886jLrbKRnuCSFizUgrI8YSKuefih\nlCwMuHK4gdHgqSgbyxMmIjJ8SUxZ20wGavZkpCsbpTh63GLCzJk4//Yxd+5NpIKzo0lTRcUoyssv\nAwDKyy+joiLYAdjUDF3lrOYrQe3zf/whIsOXAIyNQ14uLgGAq9mmQel7MQMCt9w0OCvb9MYiZvtX\nFLtPxSIyiEGYNMViQFPTRRw8OBD4UrTIXauG7q5LWcpz/sBbltYAK5cK2X4qkYG2GGmDI3t2cykS\n+YiucvSKFSsQu/KhUVRUhJ07dyZfe+mll/DGG28AABYtWoSHHnrIhmaS22IxBLsEfYWl04TSZSjz\njlhZFqO8nsAxYcttMdAGof2bidkyOZFLNIPw0NAQZFlGQ0PDuNe6u7vx85//HHv37kVOTg7q6+tx\n5513orKy0pbGEtlO9Axd0RONFNezFNAFt8XIz9g1AxpQObaSxxuSh2mWo9va2jA4OIh169Zh7dq1\naG1tTb72e7/3e/inf/on5ObmIhKJYGRkBJMmTbK1wUS247ihvWzsX0+djkWkQ0SWZTnbN7S3t+Pd\nd9/FqlWr8PHHH2PDhg04dOgQotGrSbQsy/je976HgYEBfPe73816w5GRy4hGc8W0nogonSQBJ08C\nVVV8kCLP0yxHl5aWYvbs2YhEIigtLUVBQQF6e3sxa9bYIvehoSFs3boVU6dOxeOPP655w76+i9Zb\nrRCP56O3t1/oNb1Kz4YZYeoPPdgfqULTH9d/CRiUgcHs7zU0/aET+yOVyP6Iq+zOplmO3rdvH3bt\n2gUA6OnpgSRJiMfjAMYy4L/8y79ERUUFvvvd7yI3lxmuXRIbZixdOhVLluTh6NEcTv4URfRsWs7O\ntYb9RyGimQmvXLkSW7ZsQX19PSKRCHbs2IGGhgaUlJRgdHQUv/rVr3Dp0iX853/+JwDg29/+Nm66\n6SbbGx42yg0zOjtzUVc31XO7V/lya0vRs2k5O9ca9h+FjGYQnjhxIp5++umUr82ff3Wjgf/93/8V\n3yqdJAn46CNgxozg/ztNbJiRCMSAt3av8uvWlqJn07o2O9foFpIexdnNFDa+3awj8aH/5S+P/dfN\nypUkAS0t9paHExtmNDYOoKzMe7tXZdra0tOulDxHikqEzqYdKSpJ3RnLidm5bu9CJRBnN1PY+Hbv\naDv3M1Yrq2b6upMZYCwG1NSM4s03L3qu7KvM1L30cJBResmz8Q0xG19IEqbVLR87IKC4GH2NbziS\nlQYqe7R5HTGR13g8XVFn137GyglQygxb7etuZICJ3au89Pnkp60txwWtM11C1q2mn1kbPdNlua16\n2JI9ujk5iuu0KUR8G4QTH/rHjkHoh75aUFX7eqgON9DgxYeDTOwqebpWSr2SPVrZlzpFgMrbRF7n\n23I0MPZZU1oK9PaKu6ZaWVXt64mHAa+VhykLu0qebpZSBW6PGajyNpHH+ToI20EtqGYLtmYPN3By\nSY8vlw/ZycbD4/0esGw/ZMGIgMz6JlLDIJyBWlAVeZKQkxO6/Lp8iFzilclRXDNMIeDbMWG/c3JC\nl++WD5H7PDA5ypGzh4lcxk9jlzg5oYuTx8hzdMy+5pphCgOWo13i5IQuTh4jT9FbZvZKWZzIRsyE\nXeTkkh6/LB+i4DNUZvZAWZzITgzCIeTENptEalhmJrqK5eiQ4Uxpch3LzERJzIRDRm2mdKCzY55P\n6z0sMxMBCFEQDnSQMSDTTGm1fbEDgVswEpGHhSIIBzrIGJTpoIVs64h9//By8iTXmhKRZ4UiCLu1\nWYVXA1j6TGm1dcSBeHipqgrPJCCW3Yl8J3BBOFPgc2OzCj8FMLVjCO1+eHHkIUX0CUN2shJEWXYn\n8qVABWG1wOfGWbd+2yoy0zpiOx9eHH1I8cMkIItBlFs8EvmTtyODQdkCn9ObVYgIYG6Xs+18ePHb\nQ4rdrAZRrr0l8qdArRNWO/PXDWa2ilQeNwh4Yz2vyJOjlLz0d+UFlo8P5NpbIl8KVBD22h7JRgJY\n+iYau3f/blymaEcwdIvX/q5cJyKIBuAsY6KwCVwN0K97JKeXZwFkLWe7XaoWwa9/V7bxw9g1EQkV\nuCCshxcCWHob0seQq6tHVcdj/TTzWskL/U5E5CWBKkfr4YW9k9XakKk8m6kEnWlSk9dL1V7odyIi\nrwlFJqzMwNyelStJwKuvRjO2QW95Vm3rSS9nmU72u9f7gogoIfCZcHoG1th40bVZucq2TJggY3g4\nYqoN6VkzcPW6lZXAgQP2DisqZ3HrvY9Ts6HT/75PnLDlNkREQgQ+CKdnYGfO5Lg2K1fZluHhCP7x\nHwdx990jZifCJkvQLS1Xr9vWBl3laTOBNPFzZsrKRmdDm21f+t/3yZPA9dfr/3kiIicFvhydqXQr\nalau0bJnelvMBuBs162shGaWaWViV3qQa23N0d0HevvdSvvS+7iqSv/PEhE5LfCZsF3rUc1khHa1\nRXndmpqpGBzM/v1WJnYpy8plZZexadNkdHaKm2ylNmaut33j+zhfsz+IiNwS+CAM6Ns0w2j502wg\ns2sHKmWWqRV0rIzPKoPc4CBQVzcVgJhZ2iLHzL0+W5yICAhBOVqLJAFHj+ZgyRJj5U83TmYSxeqe\n0IkgV10ttg8yjZlzKRMRBVkoMmE1yswrQW9G58VtFyUJ+OgjYMYM7dnRWtminsqA6D5Iz9BFjZkT\nEXlVqIOwMvNKMJLROVH21Fsmv/pAAZSX51nKII2Md4vsAy8+2BAR2SnU5WhlSbms7DIaG507b1iL\n0TK5yM0w3NzQhPtJE1GYhDoT9mrmZaZMLnIzDB4zSETkjFAHYcCbM2nNlMkTDxTnzuVjxgxr2bxX\nH06IiIIm1OVoJxnZ2MNsmTwWA26/XcyWlWplYRH7MnNvZyKiMaHPhO0mSUBra46hTS28momKOAmJ\npykREV3FTNhGiYBTVzcVnZ3GJjp5cYKSiAlbVq6RKYNmVk1EfsYgbCOrS6C8RsQGJWavkWk/aSt7\nTBMReQHL0TZK32d5z57fobraW9mtESLK5GavoZZBm91jmojIC3QF4RUrViB25dOyqKgIO3fuTL72\ns5/9DP/2b/+GaDSKBx98EH/0R39kT0t9yKtju1aImE2uvIbezUjUlk1xKRUR+ZlmEB4aGoIsy2ho\naBj3Wm9vLxoaGvDKK69gaGgI9913H/7wD/8QEydOtKWxfuTFJVBeYXRnrkwPNEF7yCGicNEcE25r\na8Pg4CDWrVuHtWvXorW1Nfnae++9h5tuugkTJ05Efn4+SkpK0NbWZmuDM+HkHH8yOkkr02Q1L05g\nIyLSSzMTnjx5MtavX49Vq1bh448/xoYNG3Do0CFEo1FIkoT8/Pzk906dOhWSRiScNi0P0Whu1u8x\nQpKAZcvy0dY2dqB9c7OYdbJ+IEnAyZNAVVXqe47H8zN+3+zZwCefjP9+t9TUjP2dJf7uamqm2tKu\n9P4IO/ZHKvZHKvZHKrv7QzMIl5aWYvbs2YhEIigtLUVBQQF6e3sxa9YsxGIxDAwMJL93YGAgJShn\n0td30XqrFT76aCwAA2Mf5kePDggt/xo9Z9gpaqXceDwfvb39Gb9PeUavV9bnHjhwtX8HB7XPQjYq\nvT/Cjv2Riv2Riv2RSmR/qAVzzXL0vn37sGvXLgBAT08PJElCPB4HANx4441oaWnB0NAQ+vv70dnZ\niblz5wppsF5VVbDtXF/lEphFi/LQ0yPs0pbpLeWmn9Gr9f0iS/t6rsVyMhGFmWYmvHLlSmzZsgX1\n9fWIRCLYsWMHGhoaUFJSgsWLF+P+++/HfffdB1mW8fDDD2PSpElOtDsp04QdUdmrMoB1d+di2bKp\nePvtAU8EDL2HLCi/T5kJZ/p+kbtZcWcsIiJtEVmWZSdvKLrUka38KiKQLFqUh+7uq2PYBw8OoKJi\n1BMl6kwPG5nKJ4nvKyoaxZkz6u1uacnB0qVTk38+eNB8aV/ktaxgeS0V+yMV+yMV+yOVJ8rRfqNV\npjVSbo3FgAMHLqK4+Oqa1KKiUc/s0qS3lJv4vpkzs3+/kd2slP2YqU9F7K6ldV8iIr8L3I5Z2cq0\nZrLkmTOBt98eSGacmYK8MsPz6kQuPfRuLqLsx7KysUCbfjiFHRuVsMRNREETuCCc+PBvbR2f5GsF\n0GzXTHyf6CDvtvSHBj2biyj7MXEwBTC+T0VvVGL274+IyKsCV45O2Lx5MurqUkvGIkqkiSB/8OD4\nM35FnDLkJLMHIKSfd5zIhu3eOtKuEjcRkVsClwkD6hmTqBKpWoand8ayV1ipDCj7MXEtu0vwQdyL\nm4jCzdupmknZMiY716Vmy5LtpHbO7tGjOTh6VH0Sk5XMUtmPIvpU71nBXFdMREESyExYT8ZkxwQq\nrWvadU/lOPSJE2NfW7IkLzleW1Z2GW++Of6hwCuZZaaxdCD1a42NF7Mur5IkJOcB+Pm4SCIKl0AG\nYSD7pCA7JlBpXdOuSVvpJeWTJ4G+vpyUCVOdneqlZj2Tp+ye8a3nrOBly6aiuztHtW/1PHQQEXlN\nYMrRautHM33djglUWte0a9JWekm5qmrsa4nJUsBYUDI7Pm128pYRmcriyq8VF19Gd/fVwJypbzM9\ndBC5RRqW0NLTDGmYC9opu0BkwmpZptrX7ZhApXVNuyZtjS8p52NwEHjzzYtCyrNOLAvKdlZwa2sO\nBgeBbdsmJ9ciZ+rbsrLLKZmw1yfFUXBJwxJq996BjgunUV4wF02rDiM2ITbue9rPn0JF4bxxr1G4\nBCIIqwUKu2dJK2ldU+t1KyXfTCXlWAyoqbH/4SKd2fehVhbfvHlyclOQxsaBjA8UsZi4hw4iq9rP\nn0LHhdMAgI4Lp9F+/hRunnlr8nU9QZrCIxA1O7VZvlZnSRvdIlHrmmqvO1HyNcvIjG/R7yN9U5Ap\nU9TPQU48dNTUMACTuyoK56G8YOw0ufKCuagonJfyeqYgTeEViExYLcu0kvE6ufuV2ztBaWWvene+\nEv0+/LbumggAYhNiaFp1WLXcnAjSiUw4PUhTuAQiCAPqgcLs1ol6A4qImcNuBhuRDxui34dXllAR\nGRWbEEspQae/li1IU7gEJgiLojz2TyugiApgymBTVOTsMYkis1e7xtq5PzQFTbYgTeHCIKyQHlS1\nNogQHcAqKkaFBHUj2bkd2asTQdPPp1URESUEYmKWKOlB9cyZnKwTrUQfKJBtLbHeSWJGJ0e5tdWm\nFV6eyEZEZASDsILRoCo6gKnd30jQOXkSqoE82/vw037MfjutiohIDT+9FMwEVZEBTO3+RoJOVRVc\nP+7P6NIuoz/HIw2JKCg4JpxGxJim6I03jIzbuj2j2OxkNSM/l3iPic05iIj8ip9igtkxXmk0Q3ez\nvGy2VGzm5zZvnoy6Oo4LE5F/MQgLZtd4pV/Gbc2Wio3+HMeFiSgIWI4WLOy7PJkthxv9ubD3MxEF\nA4OwYG6PyYpgdQ2u2XF1Iz8XhH4mImIQtoGfd3lycs9sq/zcz0REAMeEARhbUmN2+Y1fcKyViMg5\nof+ENTKbOQw7NXENLhGRc0IfhI1kfmHIEv24jSURkV8FL4oYZCTzC0uW6JflUER0lTQsoaWnGdJw\nAEt0ARb6iVlGZtlyRi5ROEjDkq/O+5WGJdTuvQMdF06jvGAumlYd9kW7iZkwAGOZH7NEomBLBLSl\nryxG7d47fJFZtp8/hY4LpwEAHRdOo/38KZdbRHoxCBMRKfgxoFUUzkN5wVwAQHnBXFQUznO5RaRX\n6MvRRERKiYCWKO36IaDFJsTQtOqwr0roNIZBmIhIQU9Ac3LMWO+9YhNiuHnmrba2xS1+G6M3gkGY\niChNtoDm5CSo9Hs1fuMNnOnvCmQwUhP0SWccEyYiMsDJMeP0ey17ZbGvJoyJ4McxeiMYhImIDHBy\nEpTyXsWxYnT3dwEIZjBSE/RJZyxHExEZ4OQkKOW9ivJLULd/ua8mjIkQ9ElnDMJERAY5OQlKeS89\nwSiIk5iCPOmMQZiIyCe0glHQJzEFEceEiYgCIuiTmIJIVxD+7LPPsGjRInR2dqZ8/ec//zlWrFiB\ne+65B/+uaknbAAAUdklEQVT6r/9qSwOJiLzOrcMT0u8b9ElMRvnhUAvNcvTw8DC2bduGyZMnj3vt\ne9/7Hl5//XXk5eVh+fLlWL58Oa699lpbGkpE5EVaJWC7xmjV7hvkSUxG+KU0r5kJ7969G/feey9m\nzJgx7rWKigr09/fj0qVLkGUZkUjElkYSEXlVthKwnYdBqN03MW4sIuD4IZNU45fSfNZMuLGxEYWF\nhViwYAFeeOGFca+Xl5fjnnvuwZQpU7BkyRJcc801mjecNi0P0Wiu+RZnEI/nC72e37E/UrE/UrE/\nUlntj5prb0Pl9Eq0fdqGyumVqJl7G2ITxwLgR2c+SAkE50a7UBq/3XKbte5rRaI/pEsSFr741eT1\nmzc0C7m+U0T1j93/XiKyLMtqL65ZswaRSASRSASnTp3Cddddh+effx7xeBxtbW3YuHEj9u7di7y8\nPGzatAlLlizB0qVLs96wt7df6BuIx/OFX9PP2B+p2B+p2B+pRPWHWsnZ7pKo6FK3sj9aepqx9JXF\nydcO3vOW6WVCbi2bMnNf5c+UfmGWsH8vasE8ayb88ssvJ////vvvx/bt2xGPxwEA+fn5mDx5MiZN\nmoTc3FwUFhbi888/F9JYCg5JAtrbc1BR4ewZzIn71tQ4d08KL7WlQ3aP0dq5flbUaVJOj82mB14j\n/ZPe1hMPttjWzgTD64Rfe+01XLx4EatXr8bq1atx3333YcKECSgpKcGKFSvsaCP5lCQBtbV56OjI\nRXn5ZTQ1XXQkECvvW1kJHDgARx8AiJTc3mjCbBYq6gEi09isyP5Qvj8AlgJ+eltPnjuJ6yd9SVhb\nM9EdhBsaGgAAZWVlya/V19ejvr5efKsoENrbc9DRMTb+39GRi/b2HNx886ij921rg2P3JfIaq1mo\niAcIO89nTn9/uxd931LAT29r1YwqDP5WdcRWCO6YRbapqBhFefnlZCZcUeFMIFTet7ISjt2XyGvs\nzkL1sLMkn/7+AFgK+OPaOjGGQdg7h4JBmGwTiwFNTRcdHxNW3remZioGB525L5HX2JmFGmFXST79\n/VXPmG854Ds9fMAgTLaKxeBKKThx31gMDMJkipcOQnB7XNer1N6fnw57YBAmIkrj5IxerQDrhXFd\nL/P7++MBDmQ7SQJaWnIg+W/THQopp3Zb0rOjlpW2+GnHKz+1VSQGYbJVYrnQ0qVTUVubx0BMvmDn\nQQjKYKMnwJptix1bZtoRKKVhCUfPHsGSvQtt2d7T61iOJlu5tUyJyAq7xlLTS8uN33hDc+KU2baI\nnhltpCyudwxbec0Et2Zxu4VBmGzl1jIlIqvsGGtMD4xn+rt0BVgzbRE9M1pvUDcSrJXXTCi7dg4G\nRwYhDUuBm0iWCcvRZKvEcqGDBwcc2zGLyEmJEm3PxR7NUm16abkov8TWLS2bVh3GwXveEjKxTG9Z\nXE+JPdFnRfklyWuWXTsHLy/fC0SAulfvCk1Zmpkw2c6tZUpEdlNmfRNyJmJ49FLW7E9ZWi7KL0Hd\n/uW2zsAWmc3rLYtrZeCZSvJn+rtQUTgP7edPofPChwDcKUu7sSyNmTARkUnKrG949BIA7RnMicB4\npr9Ld8bolYww21nFibYCyJqBZyrJJ65p54Q4LXae/ZwNgzARkUnKoDEhZyIA/cFDK+C4FRTMSG8r\nANVgrXzfxbFiFOWXJF8TXUI3wqllaekYhIko9LQyTrXXlUHjxNqThoKHVsCxIyjYlVkbaWtsQgyN\n33gDxfkl6Ja6Ubd/eUp7smXbdsg0Pu1kFs4xYSIKNa3ZvFqvK8ddZ+bNNHTvbGO2omc327kLmNG2\nnunvQnd/FwBxy6fMjOX2XOzBslcWo7u/a9z4tFMPAQzCRBRqWktv3DqJSPRaZTvfh9G2inzAUD5c\nlF07B3vueAbVM+ZrtkEalrBs31fRLXUDSB2fdhLL0UQUalpjs05PFlKWjEWWZu1+H0baKmrsVxqW\n8OqHjcmHi87ffqh7eVP7+VPJAAwAxfklrpwyxUyYHCVJcPxoQ6JstLI4J08i0lMyDsqJSmaWTynf\nO4Bxy8MS9GT6ymy8OFaMA/e85UqfMBMmx3AfafIqrSzOqclCredOZJ3gZHXGdGIZUPv5U56aba1n\nwlj6e1f21fDoJexa8DTKCuYA0JfpK7Pxt+uPGx7PF4VBmByTaR9pIhojDUvYdHhj8s9lBXPGBRKr\nM6ZFL3sSMdtab5vS3zuAlPL6n1bW481VRwzPUHdyJnYm/BQkxyT2kQbAfaSJ0rSfP4XO336Y/POe\nRc+MCw5Wx3VFLnsSFdD1tin9vVfPmD9uXNkLQdUojgmTYxL7SHNMmGi89BnD1TPmj/seq+O6yntY\nPShB1GxrvTOl1d67309bisiyLDt5w97efqHXi8fzhV/Tz9gfqdgfqdgfqaZcG8HR079SDWhO7yXs\nxP2kYQmt505g09sb0Xnhw5QJYEZ+P0SuOzbyvp38OxH57yUez8/4dWbCRBRK0rCEhS9+FW2ftpna\npMMOdhyfmOkeU6JTLB+UIHK2td737cbfid04JkxEodR+/hTaPm0DkHksUvS2kV46jEHPvtXKtmbb\nttPJMVi39ne2EzNhIgqlisJ5qJxemcyE1TbpEL2rkxcyuGxZbKajBu0+clEv0Vt5egGDMBGFUmxC\nDM0bmlXHhEWWW93a+jIbtRJwelv/45MmoW23MqbrtQ1HRGA5mohCKzbR+CYdZsrKbp6Ta1R6W++c\nXavZdr19ImJZkx+XIWXDTJiISCezZWU/ZXCZ2pqt7Ub6xIsVAbcxEyaiUJKGJRw/c9xQNqZnYpBX\nJjFZkd7WbG03MlnKTxUBpzATJqLQMZvRak0M8toELCfonSyVGAt248zeTO3wSkWCQZiIQketLKr1\nAa1Vmg1SuVVvsNJTanfz4UTt5CWvPCSxHE1EoZOpLKp30lC20qzyusptId1gZV2y0QlUWqV2t9b3\nZjt5ySvrjBmEiSh0EtnbsfXHktmQiECRuG7j3a8DEeg+YF40q7OQRQdNt8aCtU5e8sKYNIOwyyQJ\naGnJ4dm6RA6LTYjh9qLbk9mbqEChti2kk6wGUdFBU3l2r9ESsJWMXs/JS27jmLCLEofcd3Tkorz8\nMpqaLvJkISKXiFxG5PbOTlbvr2dZktF+MrMvttWxZD+cvMQg7KJMh9zffDPP2CVyi6gDFNI//AGg\npafZsRm5Ih4o1PrCyUlWIia6OXEohhUsR7uIh9wT+ZdWmVT54W91lygz7FqXrLfULeLAijCsK2Ym\n7CIeck/kT37cJUrU+lg9pW4R2bLIdcVeWxusxCDsslgMLEET+YyRwOr2+DAgtoSsp9Rt9cFDZHuV\n1yq7dg723PEMqmfMFzLOLQLL0UREBhkpk1qZGSxKthKymbKxVqnbahlZ5BIp5bU6f/thxmVjIg6W\nMIuZMBGRQUYnPrk9OUgtG0/POE882CLkflYnhhXll6A4vwTd/V2WqwfK957QceE0Ws+dwJToFFQU\nznN1yIBBmIjIBLcDqxFqQTE9+DSfbcbQgCykJGu2f6RhCXX7l6O7vwvFsWI0fuMNS21JvPfWcyew\n6e2N6LzwIcqunZP8//KCuWj8xhuuDRnoCsKfffYZ6urq8M///M8oKytLfv29997Drl27IMsy4vE4\n9uzZg0mTJtnWWCIiMidTUFRmiWXXzsEDbzyA059lHod1asxU+WDQLXXjTH8XZubNTPkeo22JTYih\n5osL8eaqI2g/fwqDI4Ooe/UuAGMPH2f6u1w7alIzCA8PD2Pbtm2YPHlyytdlWcZjjz2GZ599FrNn\nz8bevXtx9uxZXH/99bY1logoiNyaFKTMkNMDk7Ik6+TaYDtPqko8iEjD0rh7uFXZ0JyYtXv3btx7\n772YMWNGytd//etfo6CgAC+99BL+7M/+DBcuXGAAJiIyyOqkIKvrcRPBp3rGfFROrwQwfjKVkwcw\naE1kE7nHtxe2r8yaCTc2NqKwsBALFizACy+8kPJaX18f3nnnHWzbtg0lJSV44IEHcMMNN+ArX/lK\n1htOm5aHaDTXessV4vF8odfzO/ZHKvZHKvZHKrf746MzH6QElXOjXSiN367rZ6VLEha++FW0fdqG\nyumVaN7QjNhEcwEljnw0b2jGyXMnUTWjKuU6Ndfehsrplcn71My9zfR99Lal9AuzMr4mqi3Z7pHy\nfTb/fkRkWZbVXlyzZg0ikQgikQhOnTqF6667Ds8//zzi8Tg6OzuxceNGvPbaawCAl156CcPDw9iw\nYUPWG/b29gt9A/F4vvBr+hn7IxX7IxX7I5WyP9wqCVspr7b0NGPpK4uTfz54z1uWSqrZfj+8tOGF\nU20R+e9FLZhnzYRffvnl5P/ff//92L59O+LxOACguLgYAwMD+OSTTzB79mz8z//8D1auXCmksURE\nTnLz0Hkry3mc3AjES7PBvdQWqwwvUXrttddw8eJFrF69Gk899RS+853vQJZl3HTTTbjjjjtsaCIR\nkb3c3lrSbFARefKTXbyUQXuR7iDc0NAAAClLlL7yla9g37594ltFROQgL2wtaZaXssL0gOtmhcEv\nuFkHEYWeHzJKPdzMOjMFXLcrDH7AvaOJiKC9H7KIo/ns5Ob+x0Dmkn4YjiK0ipkwEZEGP5RVRZxc\n9NGZDzAjp8TUe8tU0g9KhcFODMJERBr8UFa1Mq4t4iFDLeB6aczaixiEiYg0+GHilpWsU9RDBgOu\ncQzCREQa/FJWNRsE/fCQEVQMwkREOgQxy1POpm5adRjnRrtMjwmTOQzCREQhlGkc+Pai27mtqcO4\nRImIKIScPBmJ1DEIExGFENfwegPL0UREIeSXyWZBx0yYiCiktHYJ8yuv726mxEyYiIgCww+7mykx\nEyYiosDw24QzBmEiIgoMv004YzmaiIgCw28TzhiEiYgoUPy0uxnL0URERC5hECYi8hg/LbEhaxiE\niYg8JLHEZukri1G7946MgZhBOjgYhImIPERriY2eIO0EPgiIwSBMROQhWktsvLAO1isPAkHA2dFE\nRB6itcQmEaQTO0K5sQ4204OAX2Yjew2DMBGRx2RbYuOFdbBeeBAICgZhIiKfcXsdrBceBIKCQZiI\niAxz+0EgKDgxi4iIyCUMwkRERC5hECYiInIJgzAREZFLGISJiIhcwiBMRETkEgZhIiIilzAIExER\nuYRBmIiIyCUMwkRERC6JyLIsu90IIiKiMGImTERE5BIGYSIiIpcwCBMREbmEQZiIiMglDMJEREQu\nYRAmIiJySdTtBpg1OjqK7du3o729HRMnTsSTTz6J2bNnu90sx61YsQKxWAwAUFRUhNWrV+Opp55C\nbm4uampq8NBDD7ncQme8++67+Id/+Ac0NDTgk08+waOPPopIJILy8nI8/vjjyMnJwQ9/+EMcPnwY\n0WgUW7duxY033uh2s22j7I8PPvgA3/rWt3DdddcBAOrr67Fs2bJQ9Mfw8DC2bt2Ks2fP4tKlS3jw\nwQcxZ86c0P5+ZOqPWbNmhfb34/Lly/i7v/s7/PrXv0Zubi527twJWZad/f2QfaqpqUnevHmzLMuy\n/M4778gPPPCAyy1y3u9+9zv57rvvTvna17/+dfmTTz6RR0dH5W9+85vyyZMnXWqdc1544QX5rrvu\nkletWiXLsix/61vfko8dOybLsiw/9thj8r//+7/L77//vnz//ffLo6Oj8tmzZ+W6ujo3m2yr9P74\n2c9+Jv/0pz9N+Z6w9Me+ffvkJ598UpZlWe7r65MXLVoU6t+PTP0R5t+PN998U3700UdlWZblY8eO\nyQ888IDjvx++LUe3tLRgwYIFAIDq6mq8//77LrfIeW1tbRgcHMS6deuwdu1aNDc349KlSygpKUEk\nEkFNTQ1++ctfut1M25WUlOAHP/hB8s8nT57EbbfdBgBYuHAhfvnLX6KlpQU1NTWIRCL4whe+gMuX\nL+P8+fNuNdlW6f3x/vvv4/Dhw1izZg22bt0KSZJC0x9f+9rX8Nd//dcAAFmWkZubG+rfj0z9Eebf\njzvvvBNPPPEEAOA3v/kNpk+f7vjvh2+DsCRJyTIsAOTm5mJkZMTFFjlv8uTJWL9+PX7605/i7//+\n77FlyxZMmTIl+frUqVPR39/vYgudUVtbi2j06siKLMuIRCIArvZB+u9LkPsmvT9uvPFGPPLII3j5\n5ZdRXFyM5557LjT9MXXqVMRiMUiShL/6q7/Cxo0bQ/37kak/wvz7AQDRaBSbN2/GE088gdraWsd/\nP3wbhGOxGAYGBpJ/Hh0dTfngCYPS0lJ8/etfRyQSQWlpKfLz83HhwoXk6wMDA7jmmmtcbKE7cnKu\n/lon+iD992VgYAD5+fluNM9xS5YswQ033JD8/w8++CBU/fF///d/WLt2Le6++278yZ/8Seh/P9L7\nI+y/HwCwe/duNDU14bHHHsPQ0FDy6078fvg2CM+fPx9HjhwBALS2tmLu3Lkut8h5+/btw65duwAA\nPT09GBwcRF5eHrq6uiDLMo4ePYpbbrnF5VY670tf+hKOHz8OADhy5AhuueUWzJ8/H0ePHsXo6Ch+\n85vfYHR0FIWFhS631Bnr16/He++9BwD47//+b1RVVYWmPz799FOsW7cOmzZtwsqVKwGE+/cjU3+E\n+fdj//79+MlPfgIAmDJlCiKRCG644QZHfz98mzouWbIE//Vf/4V7770Xsixjx44dbjfJcStXrsSW\nLVtQX1+PSCSCHTt2ICcnB3/zN3+Dy5cvo6amBn/wB3/gdjMdt3nzZjz22GP4/ve/j+uvvx61tbXI\nzc3FLbfcgtWrV2N0dBTbtm1zu5mO2b59O5544glMmDAB06dPxxNPPIFYLBaK/vjxj3+Mzz//HD/6\n0Y/wox/9CADwt3/7t3jyySdD+fuRqT8effRR7NixI5S/H3/8x3+MLVu2YM2aNRgZGcHWrVtRVlbm\n6OcHT1EiIiJyiW/L0URERH7HIExEROQSBmEiIiKXMAgTERG5hEGYiIjIJQzCRERELmEQJiIicgmD\nMBERkUv+HxaoBxF11zRcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x173ec214a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "centers = [5, 5.3, 4.7]\n",
    "std1 = 0.1\n",
    "colors = 'brg'\n",
    "\n",
    "data1 = []\n",
    "for ii in range(3):\n",
    "    data1.append(stats.norm(centers[ii], std1).rvs(100))\n",
    "    plot(arange(len(data1[ii]))+ii*len(data1[0]), data1[ii], '.', color=colors[ii])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAFJCAYAAADaPycGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtwVdXZP/BvkhMUcqKBmjDtQHy55II4lnJRO4NKX0pj\nrJdfUyiIRR1tR6mdaksVpIq0IOr08trxhrba6bTaVoTWXqSpxSJFhWaC2CEmh4AdAWUgNEGzQzSX\ns39/hBPOOZzLvqy99lp7fz8znYpAzjor8Tx7PetZzyowTdMEERERSVPo9wCIiIjChsGXiIhIMgZf\nIiIiyRh8iYiIJGPwJSIikozBl4iISLKIrBfq6OgW+vVGjx6Frq4TQr+mzjgfqTgfqTgfp3AuUnE+\nUomcj/Ly0qy/Z3nl+9Zbb2HJkiUAgHfffRfXXnstFi9ejPvuuw/xeNz9KG2KRIqkv6bKOB+pOB+p\nOB+ncC5ScT5SyZoPS8H3Zz/7Ge655x58/PHHAIAHHngAd9xxB5577jmYpoktW7Z4OkgiIqIgsRR8\nKysr8cgjjwz/uqWlBRdeeCEA4NJLL8Xrr7/uzeiIiIgCyNKeb11dHQ4dOjT8a9M0UVBQAAAoKSlB\nd3f+/dzRo0cJX87nyqeHEecjFecjFefjFM5FKs5HKhnz4ajgqrDw1IK5p6cHZ511Vt6/I3pDv7y8\nVHgRl844H6k4H6k4H6dwLlJxPlKJnA8hBVfJzjvvPOzcuRMAsG3bNsycOdPZyIiIiELIUfBdvnw5\nHnnkESxcuBD9/f2oq6sTPS4iIqLAspx2HjduHJ5//nkAwIQJE/DrX//as0EREREFGTtcERERScbg\nS0REJBmDLxERAKPfQPORJhj9ht9DoRCQ1tuZiEhVRr+Bug1z0H58L6rKqtG4YCuixVG/h0UBxpUv\nEYVerLMV7cf3AgDaj+9FrLPV5xFR0DH4ElHo1YyZgqqyagBAVVk1asZM8XlEFHRMOxNR6EWLo2hc\nsBWxzlbUjJnClDN5jitfIl0YBiLNTYDBgiAvRIujmDF2FgMvScHgS6QDw8DoujkYXT8Xo+vmMAAT\naY7Bl0gDkVgrIu1DBUGR9r2IxFgQRKQzBl8iDQzUTMFA1VBB0EBVNQZqWBBEpDMWXBHpIBpFV+NW\nRGKtQ4E3yn1JIp0x+BLpIhrFwIxZfo+CiARg2pmIiEgyBl8KPh7RISLFMPhSsPGIDhEpiMGXAo1H\ndIhIRQy+FGg8okNEKmK1MwUbj+gQkYK48vWJYQDNzYXcgpQhcUSHgdc9Fq8RCcHg6wPDAOrqRqG+\nvgR1daP4OZYBH04UxOI1ImFCFXxV+UCPxQrR3l4EAGhvL0IsFqpvQ158OFETi9fUZvQbaD7SBKOf\n/8HoIDSf+ip9oNfUxFFVNQgAqKoaRE1N3L/BKIgPJ2pi8Zq6jH4DdRvmoH7jXNRtmMMArIHQfKqp\n9IEejQKNjSeweXMPGhtPcCsyDR9OFHWyeK1r8xZ0NW7lHrpCYp2taD8+lJVoP74XsU5mJVQXmmrn\nxAd6e3uREh/o0SgwYwaDSiaJh5NYrBA1NXF+xquE/aWVVDNmCqrKqtF+fC+qyqpRM4ZZCdWFJvjy\nA10vfDghsi5aHEXjgq2IdbaiZswURIv5Aae60ARfgB/oRBRc0eIoZoxlVkIXodnzJSIiUgWDLxER\nkWQMvkSZsJMTEXmIwZcoHTs58eGDyGMMvkRpQt/JSaeHDz4kkKYYfInShL2TkzYPHzo9JBClYfAl\nZXpeKyPknZx0efjQ5iGBKINQnfOl0yV6Xic6f7Hd5Ulh7uSkyR3IiYeESPtepR8SiDJh8A25TD2v\n2YiEtHj40OQhgSgTpp1DjpcYkNYSDwkMvKQZrnxDjj2vPWAYXI0RUU4MvsSe1yKdrMBN7EOGsWCL\niPJj2plIIFbgEpEVjoJvf38/li1bhkWLFmHx4sXYv3+/6HERaUmXYzpE5C9HwffVV1/FwMAAfvvb\n3+K2227Dww8/LHpcRHoK+RlhIrLG0Z7vhAkTMDg4iHg8DsMwEIlw65homA7HdIjIVwWmaZp2/9Lh\nw4fxjW98AydOnEBXVxfWr1+P6dOn5/w7AwODiESKHA+UiIgoKBwF3wceeAAjRozAsmXLcPjwYdxw\nww3405/+hDPOOCPr3+no6HY10HTl5aXCv6bOOB+phM5HAI4OuZqPALz/ZEH6b8XoNxDrbEXNmCmI\nFjv73lidDxGvpbLE+5tdfSF6P7AdFjMqLy/N+nuO8sVnnXUWiouLAQBnn302BgYGMDg46Gx0RCoL\n+9Ehnd9/wB4a0hn9Buo2zEH78b2oKqtG44KtngXF5NeadPZk/HDOw5hWMT0wQTj5/dWeU4uXvvSK\n5+/NUcHVjTfeiJaWFixevBg33HADvv3tb2PUqFGix0bku4xHh0J0jZ22R6dCcONRrLMV7ceHvjft\nx/ci1und9yb5tfZ/sA8NL16Jug1zYPQHY16T31/bsTZP5zLB0cq3pKQEP/3pT0WPhTRlGAhsh6zT\nmvePq9R3JeiArpcXnPbQsHsXMHLk0PhzpAJ1UjNmCqrKqodXvjVjvPveJL9WQiLgzxirf3Fh8vur\nPafW07lMcLTn6wT3fL3l13y4uRXJy6Dt1Z5vJNaK0fVzh3+ra/MWLSqbQ7fnm5wunzQZABDZv2/o\nQWJXMzp6pXzseU72nu/uo7tw56t3YP/xfZ6numWTvefL4BsQfs1Hc3Mh6utLhn+9eXOPpVaVXl9l\n6Nl8aLoHGsr/Xk4+NKC3F6Mbrjz173fsQMfE8/wbl2Ls/mwEvfBK5H8ruYIv20tqyDCGgp4K21hO\nb0XKdJWhFkQ20QjR3rEwdubs5HnrgWnTU7qOYepUjwcZbNHiKGaMnRXIwCsTu2NoxusVoxMPPfQR\nAGDaNOvp40TQTrwPK0Fbmb1lEU00NF1B+yKR9h5XidENX7Q/Z2n3/pZHo0BvyLIAARKUlTeDr2Yy\nrRj9upEo04OAVXavMlTxocONTFXEOuwdC5dvPzn5IWX8eEQOHgTgYM7YdSwQZB6v8pomuT5KcJrm\n9YLb1HHiKkMrQVTbNHUWvIABlo4DpTykHDyIwfGVAEI8ZyEn83iV17jy1YzdFaOXnKSOdXgtKdJS\noVov4x2ysvpPP+rUtekviBw6ENo5CzuZx6u8xuCrocSK0W8yHwRUeugQJuSpUEtniDM8pAyMHSt/\nsKSEaHEUjQu2cs+XSOaDgCoPHWRDrj1dq6v/kD+kUKpEtbXuGHyJyBtWKroZWCmk9K5aISJladsX\nmkgCBl+ioPK5iQcruomyY9qZKIjSUr7Y1Sx/DKzo1lZQGlmojCtfogBKT/mipcWfgST2dBl4tZFo\nZFG/cW6grg1UDYMvUQClp3zZz5isClIjC5Ux+BLpwO7+rcgLIChUEo0sAGjfyEJl3PNVmDIXCZD3\ncp2HdXoJA4/xkAO6N7LQZb+aK19FJS4SqK8vQV3dKN46F2R5ehzzyA7Jpuu1gXb3q41+A81HmnzZ\n12bwVVTQLhKg7PIFVx7Z0Y+fH+phZme/2u/CMn6iZ6DCZfUq3V5E3sobXLl/q5V8H+oyA3PYHgLs\n7Ff7XVjGPd80qtwbG8iLBCgzK+dhuX+rjUwf6olexDLuo03seY4rrUTDH76I9uN7MensyfjhnIcx\nrWK6dqnkXNL3d+3sV/t9QxKDbxqVLqvnRQIhomtwzVUoFlK5PtRzBWYRkoP7+Oh4HDQOAgD2f7AP\nDS9eqf0F9MmyPchYvXjB78Iypp3TMN0rjgrpe0t8bsOorTyFYmGV+FDf/OUtpwU6r4/xJAf3g8ZB\njC+tTPn9IJ3bFZE29rOwjME3TSLdu3lzj28p5yDQplqbAcQxVmFnl+1DPVdgFiE9uL/05S3YdM2f\nMals8vC/C8q5Xd3PIzPtnEEi3ZtYuXHP1T6V0ve5ZAogWqZ/fZAoFEucP2YVtjVe3kebKZU6dtRY\nvLxgmxZnX+3wO23sVuhXvtlSo9qs3BSlS/o+kMd4ZKXRWYWtpEyrbl3P7eaj8/sK9co3V2WzLis3\nVWlTrR20m3ecdsNyStdCMSKfhXrlm6uRhS4rN5Ul0vdefPYLLeYK0M073Icl0kOog2+uAMvCK3Vx\nSyC7QKbRiQIo1GnnfKlRnrNVE7cEcghaGp1SyLw0QJcLCnQV6pUv4C41KvscqzbnZj0mdEsgiGd8\nA5RGp1Nk9iL2u+9xGIQ++OajSjU0U62nCNsS4Blf0ojXvYiT+0D73fc4DBh8c8gV8GTfOsRbjlKJ\nKOZicZIFQcwMKCzXRQheNpVIX+mOK60U/lphu+Qhn1Dv+eaTa28xkfpMHFPyuhpa9utp72TQyLXv\nySYRecg+tqQwGfuf+S5d8LKpRPpK91D3AaGvJeNCCd2Ee/mUh0rV0Ky+tsEwgFmz8qeT2SQiJ2YG\nhsja/7SS6vWqqUSmVbXI12Ia+3Rc+eagWjU0q6+ticRagba2oX/O1zKSTSKyYmZgSHrg2H10F0ZG\nRgpfffp5xZ3XrRr9vr5PRQWmaZoyXqijo1vo1ysvLxX+Nb1kGPC025Nu8+Epw0D5Ff8LtLWFPl2a\n4PjnI4BXBtqdi+SU6aSzJwMFwP7j+zxJn/pxvEfWZ4cuR5dEzkd5eWnW3+PKV4JcbSzJA9Eo0NSE\nru3/ClTQkCYt4IY9M5C8Kuwd6EXDi1cC8OY+Xi8vXfBbkN+bE9zzlYCVyj5wc9Y1zBW+PH6VUSJw\nTKuYrvU1dqQOrnwlYKWyRkJe4csrFnPT/Ro7UofjJdiTTz6JhQsXoqGhARs2bBA5psBhpbI+wl7h\ny97Q+dmpAubZVsrG0cp3586dePPNN/Gb3/wGvb29eOaZZ0SPK3BYqayH0Ff4sje0MDzbSrk4Wvlu\n374d1dXVuO2223Drrbdizpw5godFOghkr+l8Z3/d7AfrspfsZL9cl/cmkd9nW7nqVpujo0b33HMP\n3n//faxfvx6HDh3C0qVL8de//hUFBQVZ/87AwCAikSJXgxXBMICWFmDqVD7Uu3GyjwXa2oDaWqCp\nKQTz6eZNB3nCvHxvkv6DNfoMtBxtwdSKqYiOEPM6Rp+BWT+bhbZjbag9pxZNX28S9rVVfm2yxlHa\nuaysDBMnTsSIESMwceJEnHHGGejs7MQnPvGJrH+nq+uE40Fm4uQsVpCP/Mg+59vcXIi2thIAQ5+5\n27f3KJVW92I+Is1NGH2yeQfa2oaOMlksRnLzd0Xw8ufDs/fmUfFb+lx4mR5+6UuvDBdn9X5gohfi\nvweZzs82H2lC27Gh70nbsTZs3/uvrMd82CMglaxzvo7SzjNmzMA///lPmKaJI0eOoLe3F2VlZY4H\nKIvMIz92UrI6pm/dXuvn63t2mCJ1U4wU5EImr96brOI3L9PDXrWDTMjW+tLLSxhksps61ynV7mjl\n+7nPfQ5NTU2YP38+TNPEqlWrUFTkf0o5H1lHfuyssHVdjedrvZmL6Pdsq3uYm9WUm2KkIBcyefTe\nZBW/6dz6MNODQyLY634kym5GQrcCN8fnfO+66y6R45DCTcCwI9dtSG7+rGqcVnCLfM92A7nrc6xu\nOj7Z/bs6tXb0ohOWpAcWnQNVrgcH3TtKZXuwEPXn/Ra6Vksi7oHNx05K1m361gmrKV+vUsMi37Pd\nrQRt0r9B7DTlJN3vplOZDV6nh72SeHDY/OUtyq/07LKbOtct1c6LFTxiJxWa/medXMJgdT6srhS9\nToeLumgi2zhzzocGK8pIcxNG188d/nXX5i2uVpa+//eiUOcw3+dCMSrPh93LGERc3qB0wRXlZ2eF\nnfxnE8Gkvr4EdXWjhC94rK4UvS5OE5WBcNQ9TNRqKt9KzsXZV21W6BaFvXMYOWM3I6FTBoPBVzFe\nBz2rKV8/0uFO2QrkoppB5EsLu00b52v2oRkpDxMKNPrQpdpWl3EGGS9WUIzXFdlWi85kFadJJTD1\nma9wS8gFBUG6zs/r4ikF0tq6VNumj3PX0ma/hxRKXPkmcVJgJLIoKbEPummTt5cwWF0pyihOk0lk\n6jPfSi5oaWMhkh4mRK9QVUhr+91O0qr0cbYcbfF5RHKpsupn8D3JyV6ryP3ZI0eAyy4rQX19CRoa\nRnm22tSxoYeoMQsNiPnSwgFLGwvjURW3Cg87ulTbpo9zasVU38YiOxBma0riB6adT3Jy9lTUeVXD\nAK64YhQOHix0/bXyvY5uDT2Ejll06jNfWjhIaWNBPLsvWIEmJnbOC4uoynXqtHGOiHrS9jIfJ2l6\nt/Om0llgrnxPclJgJKooKRYrxMGDpzqEjR8f96TASWZ7TVGEj1nSuVHKzNMVqgLfWyvVtiqsvlSo\nCrabphcxbyplJ7jyPclJgZGooqTkIqvx4wfx0kverEhltdcUKdOYRZ0RJh8osEL1m0qrLz/Zbesp\nYt5U6mbG4JvESbtEpy0W07+GjMpiHSuY08cMZE5DaxuQNWj4IVzI0/FWg46fqWkZkgPhuNLKvO9V\nRA9uleaUwVcRIoI4kD8IiXodmZLH3Nx8ehq6piZ+WkAuL/dzxBYpcDyGxLDzoW5l9aXLsaVkTgJb\ntDiKmjFTLL1Xt6tW1eZU/U0/AXSs8HXC6+5YKsi0z67jXjagxvEYcs/JXmS+PVddji0luNmPtfNe\nc81bvspp1eZUj08pF8IQkBJ0DUJ2ZGonaanwTUb3I5uvocLxGMsU6B7liIRxe/GhrlJhkBVu5kDE\ne7US/FWb08CnnXW9ss/JHqaIgiod9k7TU+d597JlpHedvIab4iOZe8W6pscljduL+4BVKgyyws0c\niHivVoqxVJvTwAdfvyt8nQQzp2db3RZU6XgOOCHXXrZnZ0tFvIaT4iPJwVDG/HlB1ri9+lDX6T5e\nt3Pg9r1aDf4qzWnw8pJpHN16I4jTlLeb9LGblpBBTVvLSO/KTCHL3ivWKj2eROa4VTg36zc/50DH\ne40Dv/IF/KvwdZry9mu17vZ1lU1ZyzhbKvH8aiKoJFa+ngdDXc/m6jpuckSlVa0VBaZpmjJeSPRl\nzapdAJ0p8LhJ49oNZKLmw2kAVSplbRgoP3oAHRWV2T9wdT9fa3P8qv334ifORaqgz4fdI1Ai56O8\nvDTr74Vi5eu1bIHHzR6sX6t1p6+rTGHbyf1QtO/F6Gz7oboWECULeaMK8p9KDSuyUe1sb7JgbOr5\nLNdeadCu5ctGVJ9rt6zsh2b8M7oepXEqbO/XZ6pcYydKrqM9Kr1X1c72JmPwFUCVwOMnq4VtXjc8\nsVJkc9qfGVfpyTV3w1QLdB5d60eZqXCRgmjZgppq71W1s73JmHa2Idt+qI49k72QL2V95AhwxRUl\nOHiw0Lt94ZNFNuVHD6Ar255vWiGOp0dSFExx63p0SFdBvEgh29Ee1d6ramd7k3Hla1G+Y0Oqp5f9\nbrGZ7c5iT0SjwEUX5Q5ySdfPeXkkRcUWkroeHdKVyqsvp7Id7VHxvap6DIwrX4uUKShyQGbVdTay\n7ix2xMMjKdKPBVkRhCM4GlWrq7z6ciPT0R6d3qvfBWNc+Vqk275u8krXafMMkX2xk+dv6M7iHrU+\nM726iP1koOvavEWJlPMwBS6edyzTnrVq++ppVF19eUGH96rC3jSDr0XpBUWAujclJQfNefNGobMT\nmDTJ/oODm6CdaW4eeugjbNrUg1dfPYGxY62/H+3pEugUD2AJp6Xyd+9iARnZokIVNINvmlx7o8kF\nRSrflJQcNPfvL8J115UAADZtstdi08lqP9NqOfHvGhpKsHz5mc7eVJbXyvkApEkwUYJGFdDpe9YA\nlNtX94tKx3xUpsLeNPd8k1jdG1V9/ze5TWTC/v1FGDkyY7+JrHu6Tqq4s62WRc9X3u+VglXGKtOq\nAjp9zxpQb1/dByo3lFCNCnvT2q98RVbxWk2zyt7/tfseE0Fz06aenOlmK3u6yVXcVsaRaW68mK9M\n36uU8bW0cDVkg3YV0MmpfFX31SVTIZWaiZPVuIwVvN9701qvfEX3E7Z6sYDMc71urhecPTuOl1/O\nPk47K3ir48g2N6LnK/17NW5cPGV8u7ZNxZkiV0OqVteKGpfuFdBst+nJvcJuOVmNy17B+1X1rPXK\nV/QVeHauH5R1rtfte8w1TjsrUjvjyPSaoucr/Xt16FDq+FreFbgaUnU/1Om4su2FOy0M03BvPYh7\no3av1Rueg778c+B0vpysxmWu4P2setY6+HqRzlStWYaXKW47DxsqHrVK/l6lj2/qVAirMlaxUQbg\ncFyiHyRUfTDJwejz/5iJV6ymUpODzqyfzco5B24ClJPCJpnFUH6m6rVOO+ve1tFKAwuv36PVW4xU\nn+vTx1eK3l4xX1vJRhlwNi7RhVVaFWqd1HK0RakWiH5IDjptx9pyzoGblpFOCptkFkP5marXeuUL\nqLdStcpOAwtV3qMq48jGs/GpWtDjYFyiC6uEfD2RaWsLX2tqxVTfj5nY4UWKPHl1WXtObc45cLsS\ndVLYZOfvuJkfu6l6kQpM0zRlvJDoy5p1vwC6ubkQ9fUlw7/evLnH1fEb3edDNM5HqpT5EF085ubr\niTwSZvFrlZeX4j/vH8buo7sAANMqpvtW8Zqv2MfL4qPEa8+uvhC9H+QOA363YszGi/kR+dlRXl6a\n9fe0X/nqSsU9VNKM0xWj273w9Nd18fVE7qfb/VrLX/0OGl68MuM+poyCLCt7qekp3xf3bRI2puHV\n5Yjs37fEPABQsmWkqserrGDw9YmdYifZ/L4BiZA/sPpV6CT4dUWmwe18rVwf2rIqYK0EjuSUb3Hh\nCHz7H9+UViSmQv/jfFToVOUUg69kyYFNxT1UkZcpkEMWApxfFdhCX/dkurpr01/E7Kfb2APP9aEt\nazVlJXAk9iT/73OPoj/e5/mYkumwqvRzz9YtBl+JdAhsos9Ok02GgTNe3JQ3wPnVkUrY6yY/YDR8\nUdz+s8UUeK4PbVmrKauBI1ocxTWTG6Sv8HRZVfrdqcopVwVX//3vf9HQ0IBnnnkGkyZNyvlnWXAl\nvsgq+ajShAli5kN01zC/6PjzkVwwZBaPQEF/X+4iJBuFTkLnQ0DBVqS5CaPr5w7/umvzFmlHlKzM\nhYoFRl6NKdd8qDgPXpNVcOX4nG9/fz9WrVqFM88Ud0tN0FltX2lFepDctcv+3890Zlf187xBlpzS\nLejvw4f/9yg+vqYhe4Dzq6Wim9dNBO5xle7PTnvY8jPTRfGWh+VRwHIzJp1eMywcB9+HHnoIixYt\nwlNPPSVyPIEmMrClp4dbWoCJE6393XyrW6uNNygDFwEhvWlGzsCrsmxzkH4UaNNfEDl0wP8jSi4l\nB1sAvFmIrDEd2Lhxo/nYY4+ZpmmaX/3qV819+/bl/Tv9/QNOXoqy6O42zdpa0wSG/r+72/rf3bFj\n6O8l/rdjh3fj1F5399AEWZlgN98UJ68n6uuKfM1ccyDyB0+BH+Luj7vNV955xax+pNrEapi1j9aa\nr7zzionVGP7fjoM7Tvs7Ow7uMLs/Fvz9dUnmuFSdA9kcrXw3btyIgoICvPHGG2htbcXy5cvxxBNP\noLy8POvf6eo64fgBIZ1hAEePlqKiotvW/bR+Ez22l15CSjtFq/sUFRVAVdWplW9FxQl0dLgfj0qE\n7NvYXF1Fmpswuq1t6Bdtbeja/i9n6dmJ5wG9JtArbo+2/Ir/BdraTn8fgleQOeegohKjk1b2XRWV\nQKbvkZXsgdWvlYGIn43k5g4JbcfacPyDEyntCisKK4dfS9X7dkeeXYDpT8yQMi5V5yCZ0nu+zz77\n7PA/L1myBKtXr84ZeEU6lTIdCiDJKVOVi4W8GJvV9HB60Oe+rjV2+xar2gM6EmsFTgbE9Pchujdz\nzjmwcm2h1YcBn69ATD6Gk1BVVo1pFdOz9iV20yPZSzJ7Xas6B37Q7hxJrqMwKh+T8Wts2Y43qXjG\nWAlJzS1sH6tRtAf0QM0UoLZ26J/T3ofwI0v55iDPUSBb54gF3VrlRPIxnElnT8ama/48vIrLdvRF\n1aM7MntdqzoHftCut3OuFWTYVr7JsqVKRB9v8orolLyj1FGmVReg7wXzScpHFgylgDO9Dw+rhm3z\nspDq5PscPftCdPS6/9hzUtWs4tGdRK/rWGcrxpVW4lD3AU/Hp+IcJJOVdtYu+ALc880k2w+M7AcS\nJ+/RizE6+Q/Iz7OnXtPq3LMXDwNJQR21teh46RX/HzQUkfjZ0GE/VgZerJBDNApcdFHWrSBl06l+\njE1mD2mnHbxU2S7wq2tUILm5JtCDdHJyOhttbdLacepEh3aSQaJl8CV7ZAV9p0FUmRueZO3Ziry/\nVkV+XfqQQ/KDFWpr+WCVAfdj5XLcZIMondMOXkpVX3vdNcrtnqZK+7NZiK6gFiKpOnr07AuHjnJJ\n4tUep+ivm+g1rfJ+bJAw+JIwboJoWLpquQpMCnV1ysWzI1duHzyS09mizlDn4dU+qldfl+0k5Qll\n2jkI99Wq+h5U3nNXgZt9ZdfX+clKd3uRvlcwlW2FV/uo3J/VX+iCrw7X+uUThPcQZt0P/QRdm/5s\nOzA5DtyGAfzjHxg971J5wUtw0ZRf9xe75dU+Kvdn9Re6tHOmoiDd0p2Z3sOECe6+pspHtAIjyzli\ny5x0dTr5mmjfO/wfuzL7sDao2j0sH6/2Ubk/q7/QrXyVqax1QfR74EpaDiGrN5srypQjNifpFLyG\nKdo9zAqvLnvX9RJ5GhK64Ov23KsKe62iz+6qcs426Pw4R5zympMmO0p3K8PHdpKkJ6PfQPORJhj9\n6q0oQpd2BpxX1trtxORlKjf9PSQeCpy8ltMjQgR7Fbh+XAZw8jXLjx4YuvmHgYtCQvWOXYFc4ni1\nOrWzQpSZyjUMYNYsOH4tOytpFVb+ynBSgevH6i1XSzjyncqrM52pXhEeuOBrJ+jZDSR29lplpnJj\nscLEjXGOX8vKESHuDafStQKX1JFYndVvnIu6DXMYgAUaV1qJ8aWVANSsCA9c8LUa9JwEEjsrRJmF\nXTU18cSNcZ6+FveGU7EXtMI0aeGp+uosG9VX60a/gYY/fBEHuw9gfHQ8Nv2/vyiVcgYCFHwTq9hx\n46wFPacGFn7ZAAASc0lEQVSBxGoTCZkXGkSjQFMTPH+tIFSKCyW6AtdNwNAk2EihUUMOHc/r6rBa\nT36oOWgcxKHuAz6P6HSBKLhKL4TatOkEDh3KXXwko8hIZstEGa+lVA9mVYjqBe2mdaQmbSdtc9hO\n0qve0l70aNbxvG6m1bpqLSkTDzWJYquaMVOUu0c4EME3fRV76FD+xhkMJM6EpQezbG4ChuO/q/Il\nDS4eKLxoyOFl5ayf/ZSNfgPvHHobFYWVlt9PpsCmmvSHGgDKVT4HIu3sNB3KPsRqC1NltZv9Y0d/\n121q1uM0t6tiNg8acui6N5tL4oHi4qcvtpU+TgS2zV/eokQQyya5CYmK379ABF+Z+6skR+gqq90E\nDAd/11Vwk7Cn6rqYTfCRLh33ZvNxE5B0666l4vcvEGlngOnQoAlCD27bMu0fW00N29x7dpOalXJf\nrx8NSXINR8O92UyS9z11SB+LouL3LzDBl4IlMF233OyrellI5SK4SbvkQFQxWx5WC3F0v+s20751\n44KtOBo/YGvPV1eqff8CG3z9uqWHtwOJEYiCOJfB0/MVptPgptiq1A3VWxCKlK1KeUL5Rejo6PZ5\ndOETiD3fdH7tF4Zun1Kw9AIr3Qvi3HbAUrqJR0AuOVCxEMcrKu57hlkgg69fnZjYAcq5ID64iCga\n0vUavRQKNwCxEpBU7+ZklS5VymERyLSzX/uFgdmn9EEgC6xEpGcl7Xt6RvEGIPkKcYKWllZt3zPM\nArk08+voEY88ORfY1pUBSc86pcPlE7mOzchKSwdldU3WBXLlC9g7eiSySIpHnpxJPLjs3h3I50H5\nFOleJa0y2iMyjuMEYXWtWutGHQQ2+FqV3heaK1Z/LV9+Jr8XbmVK9ZaX+jMWzSujZZwP1aFXci5B\neHjwQ+iXGVaKpMLU5tBPLFgTw9dUb6biKoVS70a/gZ2HdtpK73rdzUn3KmSvU/NBTcmH/tMt315j\nEKtwVZXre8EHIOt8O6Kk+FV+TnsZe033KmQvHx50uL7QqcCnnfPt5+Zr5hDIKlxFZftecGvAoqR9\nXj9SvVLaTrqgcnpX5ypkL1PzKn/P3Ar0ytfqqjVXM4fAVuEqKtP3guloC9JXnYD0VK/STUGgf3pX\nZV6l5oP8PQv0ylfEqjUQbQ41x/PT+Smx6lS8uCqxQgtLL2O/iaiAVvFCBFECHXxFfWjz+JC/+ACU\nnzJHetw2BfH4iFS0OMpexhKIrIDWOSWfS6CDLz+0g4MPQHkovuq0RPFuWGRdkPdqRQn85pnuzfmJ\nLFPoSI8Tkd27lO+GRdYEea9WlECvfIlIE4aB0jvvGP7lwKTJyhVsUWaZ9naDvFcrCoMvEfkuEmtF\nZP++4V93//BhbVfwYZJrb9fqXm1YW1MGPu1MROo77ZjStOk+j4issNrdKluXqvQmGkdOHAlkN6tM\nuPIlIv8pVDDm90rM79e3w8rFE7lWx+nB+4qNc3Gw+0AoekQ7Cr79/f1YuXIl3nvvPfT19WHp0qWY\nO3eu6LFREpE3LxEpSYG7i/2+JMDv17fLyt5ursrn5OA9PjoeB7sPZPxzQeQo7fzHP/4RZWVleO65\n5/Dzn/8ca9asET0uLXnVf5j9pYnkkHV/r6qv70S+7la5Kp+T+1q/NP+VUFVIO1r5Xn755airqwMA\nmKaJoqIioYPSkZf9h9lfmkgOGff3qvz6Xsi3Ok4uzApThXSBaZqm079sGAaWLl2Kr3zlK7jqqqty\n/tmBgUFEIsEN0jt3AhdffOrXO3YAF10k5msbBjBrFtDWBtTWAk1NLARVimEALS3A1Kn8xgSA0Weg\n5WgLplZMRXSED3u+Pr8+yeE4+B4+fBi33XYbFi9ejPnz5+f986LbuZWXlyrVIs7rm3fy7fmqNh9+\nkzYfmnRl4s/HKZyLVJyPVCLno7y8NOvvOUo7Hzt2DDfddBNWrVqFz372s44HFiRet7Jke0U1KXGh\nARFpx1HB1fr16/Hhhx/i8ccfx5IlS7BkyRJ89NFHosemHbayDB/Vr9EjIjW52vO1I+hpZ79xPlJJ\nnQ+Pb+IRgT8fp3AuUnE+UslKO7PDFdnm1ZEqbWl+oQERyccOV2SL14VlRERhwJUv2ZLpzDEREdnD\nT06ypaYmjqqqQQBAVdUgampYgR0IhoFIcxO4l0AqynYxg86YdiZb7B6pYk9qDWhyVpnCSbd+11Zx\n5Uu2WT1SxZ7Uesh0VplIFTr2u7aCwZc8w/1hPfCsMqks18UMOgtE2pmpTTUl9ocTldHcH1aUQnfp\nEqWzcm2hjrQPvjz64j2nDzdet9wkgRS4S5com+Sbj4JC+zwgU5vecrtvy5abRESn0z5S8eiLt/hw\nQ3S6IB59Ibm0Tzsztekt7tsSpQrq0ReSS/vgC/C6PS/x4YYoVaajL0HbjyTvMYdIeXHfluiUoB59\nIbkCsfIlIpJF5tEXo98I3BEbGsLgGwCGAbzzDlBRwSOaRDLIOPrCveVgY9pZc4mjQBdfDLZwJAqQ\noLZVpCEMvprjUSCiYOLecrAx7aw5HgUiCqagtlWkIQy+mkscBTp6tBQVFWytSRQkQWyrSEMYfAMg\nGgUmTAA6OvweCRERWcENQiIiIskYfImIiCRj8CUiIpKMwZeIiEgyBl8iIiLJGHyJiIgkY/AlIiKS\njMGXiIhIMgZfIiIiybQNvoYBNDcX8hYfIiLSjpbB1zCAWbOA+voSXqNHRETa0TL4xmKFaGsb+mde\no0dElJ/Rb6D5SBOMfq5WVKBl1KqpiaO2duifeY0eEVFuRr+Bug1zUL9xLuo2zPE8ADPQ56flrUbR\nKNDUBGzf3oOamjiv0SMiyiHW2Yr243sBAO3H9yLW2erZVYWJQN9+fC+qyqrRuGAr7yLOQMuVLzAU\ngGfMYOAlIsqnZswUVJVVAwCqyqpRM2aKZ6+VKdDT6bRc+RIRkXXR4igaF2xFrLMVNWOmeLoSTQT6\nxMrXy0CvMwZfIqIQiBZHPUs1p7+OrECvMwZfIiISSlag15m2e75ERBROQaimdrzyjcfjWL16NWKx\nGEaMGIG1a9fi3HPPFTk2IiKiFEGppna88v373/+Ovr4+/O53v8OyZcvw4IMPihwXERHRaYJSTe04\n+DY3N+OSSy4BAEybNg179uwRNigiIqJMZB6b8pLjtLNhGIgmHbItKirCwMAAIpHMX3L06FGIRIqc\nvlxG5eWlQr+e7jgfqTgfqTgfp3AuUqk0H0afgZajLZhaMRXREaenk8tRil1Lm3P+GbdkzIfj4BuN\nRtHT0zP863g8njXwAkBX1wmnL5VReXkpOjq6hX5NnXE+UnE+UnE+TuFcpFJpPuzs50484zz0fmCi\nF2LHLnI+cgVxx2nn6dOnY9u2bQCA3bt3o7q62umXIiIiCsx+rhWOV77z5s3Da6+9hkWLFsE0Taxb\nt07kuIiIKGTC1B3LcfAtLCzED37wA5FjISKiEAtTdyx2uCIiImWEpTsWO1wREZFUQehQ5RZXvkRE\nJE1QOlS5xZVvABkG0NxcCEODh0qdxkpE7oWpojkXBt+AMQygrm4U6utLUFc3SumgptNYiUiMoHSo\ncotp54CJxQrR3j7USay9vQixWCFmzIj7PKrMdBorEYkRpormXLjyDZiamjiqqgYBAFVVg6ipUTeY\n6TRWIhInUdEc1sALcOUbONEo0Nh4ArFYIWpq4ogq/LOt01iJiERi8A2gaBTapG91GisRkShMOxMR\nEUnG4EtERCQZgy8REZFkDL5ERESSMfgSERFJxuBLREQkGYMvERGRZAy+REREkjH4EhERScbgS0RE\nJBmDLxERkWQMvkRERJIx+BIREUnG4EtERCQZgy8REZFkDL5ERESSMfgSERFJxuBLREQkGYMvEZFG\njH4DzUeaYPQbfg+FXIj4PQAiIrLG6DdQt2EO2o/vRVVZNRoXbEW0OOr3sMgBrnyJiDQR62xF+/G9\nAID243sR62z1eUTkFIMvEZEmasZMQVVZNQCgqqwaNWOm+DwicoppZyIiTUSLo2hcsBWxzlbUjJnC\nlLPGGHyJiDQSLY5ixthZfg+DXGLamYiISDIGXyKiDHikh7zEtDMRURoe6SGvceVLRJSGR3rIawy+\nRERpeKSHvMa0MxFRGh7pIa8x+BIRZcAjPeQlR8G3u7sbd955JwzDQH9/P1asWIHPfOYzosdGREQU\nSI6C7y9+8QtcfPHFuPHGG/HOO+9g2bJl+P3vfy96bERERIHkKPjeeOONGDFiBABgcHAQZ5xxhtBB\nERERBVmBaZpmrj+wYcMG/PKXv0z5d+vWrcMFF1yAjo4OfP3rX8fKlStx4YUX5nyhgYFBRCJF7kdM\nRESkubzBN5tYLIbvfOc7uOuuu3DZZZfl/fMdHd1OXiar8vJS4V9TZ5yPVJyPVJyPUzgXqTgfqUTO\nR3l5adbfc5R23rdvH26//XY8/PDDqK2tdTwwIiKiMHIUfH/84x+jr68P999/PwAgGo3iiSeeEDow\nIiKioHIUfBloiYiInGN7SSIiIskYfImIiCRzXO1MREREznDlS0REJBmDLxERkWQMvkRERJIx+BIR\nEUnG4EtERCQZgy8REZFkjjpc+Skej2P16tWIxWIYMWIE1q5di3PPPdfvYUn3pS99CdFoFAAwbtw4\nLFy4EPfffz+Kioowe/ZsfPOb3/R5hHK89dZb+NGPfoRf/epXePfdd7FixQoUFBSgqqoK9913HwoL\nC/Hoo49i69atiEQiWLlyJS644AK/h+2J5Ll4++23ccstt+B//ud/AADXXnstrrjiilDMRX9/P1au\nXIn33nsPfX19WLp0KSZPnhzan41M8/HJT34ytD8fg4ODuOeee/Cf//wHRUVFeOCBB2CapvyfD1Mz\njY2N5vLly03TNM0333zTvPXWW30ekXwfffSRec0116T8u6uvvtp89913zXg8bn7ta18zW1pafBqd\nPE899ZR55ZVXmgsWLDBN0zRvueUWc8eOHaZpmua9995r/u1vfzP37NljLlmyxIzH4+Z7771nNjQ0\n+Dlkz6TPxfPPP28+/fTTKX8mLHPxwgsvmGvXrjVN0zS7urrMyy67LNQ/G5nmI8w/Hy+//LK5YsUK\n0zRNc8eOHeatt97qy8+Hdmnn5uZmXHLJJQCAadOmYc+ePT6PSL62tjb09vbipptuwvXXX4+mpib0\n9fWhsrISBQUFmD17Nl5//XW/h+m5yspKPPLII8O/bmlpGb5X+tJLL8Xrr7+O5uZmzJ49GwUFBfjU\npz6FwcFBdHZ2+jVkz6TPxZ49e7B161Zcd911WLlyJQzDCM1cXH755bj99tsBAKZpoqioKNQ/G5nm\nI8w/H5///OexZs0aAMD777+Pc845x5efD+2Cr2EYw+lWACgqKsLAwICPI5LvzDPPxM0334ynn34a\n3//+93H33Xdj5MiRw79fUlKC7u7g389ZV1eHSOTUzolpmigoKABwag7Sf16COjfpc3HBBRfgrrvu\nwrPPPovx48fjscceC81clJSUIBqNwjAMfOtb38Idd9wR6p+NTPMR5p8PAIhEIli+fDnWrFmDuro6\nX34+tAu+0WgUPT09w7+Ox+MpHzphMGHCBFx99dUoKCjAhAkTUFpaiuPHjw//fk9PD8466ywfR+iP\nwsJTP86JOUj/eenp6UFpafYLroNi3rx5OP/884f/+e233w7VXBw+fBjXX389rrnmGlx11VWh/9lI\nn4+w/3wAwEMPPYTGxkbce++9+Pjjj4f/vayfD+2C7/Tp07Ft2zYAwO7du1FdXe3ziOR74YUX8OCD\nDwIAjhw5gt7eXowaNQoHDhyAaZrYvn07Zs6c6fMo5TvvvPOwc+dOAMC2bdswc+ZMTJ8+Hdu3b0c8\nHsf777+PeDyOMWPG+DxS7918883497//DQB44403MHXq1NDMxbFjx3DTTTfhzjvvxPz58wGE+2cj\n03yE+efjD3/4A5588kkAwMiRI1FQUIDzzz9f+s+HdkvGefPm4bXXXsOiRYtgmibWrVvn95Ckmz9/\nPu6++25ce+21KCgowLp161BYWIjvfve7GBwcxOzZs/HpT3/a72FKt3z5ctx77734yU9+gokTJ6Ku\nrg5FRUWYOXMmFi5ciHg8jlWrVvk9TClWr16NNWvWoLi4GOeccw7WrFmDaDQairlYv349PvzwQzz+\n+ON4/PHHAQDf+973sHbt2lD+bGSajxUrVmDdunWh/Pn4whe+gLvvvhvXXXcdBgYGsHLlSkyaNEn6\nZwdvNSIiIpJMu7QzERGR7hh8iYiIJGPwJSIikozBl4iISDIGXyIiIskYfImIiCRj8CUiIpKMwZeI\niEiy/w/ctMQjMIJuRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x173ec2146a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "std2 = 2\n",
    "data2 = []\n",
    "for ii in range(3):\n",
    "    data2.append(stats.norm(centers[ii], std2).rvs(100))\n",
    "    plot(arange(len(data1[ii]))+ii*len(data2[0]), data2[ii], '.', color=colors[ii])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Note:** In both cases the means have the same difference, but the variance is much larger in data2!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ANOVA with Sample Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Get and sort sample data\n",
    "\n",
    "*Twenty-two patients undergoing cardiac bypass surgery were randomized to one of three ventilation groups:*\n",
    "    \n",
    "  - *Group I: Patients received 50% nitrous oxide and 50% oxygen mixture continuously for 24 h.*\n",
    "  - *Group II: Patients received a 50% nitrous oxide and 50% oxygen mixture only dirng the operation.*\n",
    "  - *Group III: Patients received no nitrous oxide but received 35-50% oxygen for 24 h.*\n",
    "    \n",
    "*The data show red cell folate levels for the three groups after 24h' ventilation.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Get the data\n",
    "inFile = 'altman_910.txt'\n",
    "url_base = 'https://raw.githubusercontent.com/thomas-haslwanter/statsintro_python/master/ipynb/Data/data_altman/'\n",
    "\n",
    "url = url_base + inFile\n",
    "data = genfromtxt(urllib.request.urlopen(url), delimiter=',')\n",
    "\n",
    "# Sort them into groups, according to column 1\n",
    "group1 = data[data[:,1]==1,0]\n",
    "group2 = data[data[:,1]==2,0]\n",
    "group3 = data[data[:,1]==3,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Levene-test\n",
    "A Levene-test and/or a normality test should be made before applying a oneway ANOVA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: the p-value of the Levene test is <0.05: p=0.045846812634186246\n"
     ]
    }
   ],
   "source": [
    "# check if the variances are equal with the \"Levene\"-test\n",
    "(W,p) = stats.levene(group1, group2, group3)\n",
    "if p<0.05:\n",
    "    print('Warning: the p-value of the Levene test is <0.05: p={0}'.format(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### One-way ANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The results from the one-way ANOVA, with the data from Altman 910: F=3.7, p=0.04359\n",
      "One of the groups is significantly different.\n"
     ]
    }
   ],
   "source": [
    "F_statistic, pVal = stats.f_oneway(group1, group2, group3)\n",
    "\n",
    "print('The results from the one-way ANOVA, with the data from Altman 910: F={0:.1f}, p={1:.5f}'.format(F_statistic, pVal))\n",
    "if pVal < 0.05:\n",
    "    print('One of the groups is significantly different.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Elegant alternative implementation, with pandas & statsmodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                df        sum_sq      mean_sq         F    PR(>F)\n",
      "C(treatment)   2.0  15515.766414  7757.883207  3.711336  0.043589\n",
      "Residual      19.0  39716.097222  2090.320906       NaN       NaN\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Programs\\WinPython-64bit-3.6.0.1Qt5\\python-3.6.0.amd64\\lib\\site-packages\\scipy\\stats\\_distn_infrastructure.py:875: RuntimeWarning: invalid value encountered in greater\n",
      "  return (self.a < x) & (x < self.b)\n",
      "C:\\Programs\\WinPython-64bit-3.6.0.1Qt5\\python-3.6.0.amd64\\lib\\site-packages\\scipy\\stats\\_distn_infrastructure.py:875: RuntimeWarning: invalid value encountered in less\n",
      "  return (self.a < x) & (x < self.b)\n",
      "C:\\Programs\\WinPython-64bit-3.6.0.1Qt5\\python-3.6.0.amd64\\lib\\site-packages\\scipy\\stats\\_distn_infrastructure.py:1814: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= self.a)\n"
     ]
    }
   ],
   "source": [
    "df = pd.DataFrame(data, columns=['value', 'treatment'])    \n",
    "\n",
    "# the \"C\" indicates categorical data\n",
    "model = ols('value ~ C(treatment)', df).fit()\n",
    "\n",
    "print(anova_lm(model))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
